      subroutine ornl_dfir1d(	in_put, incinp, i0_inp, n_inp,
     $				firfil, incfir, i0_fir, n_fir,
     $				output, incout, i0_out, n_out,
     $				alpha, beta)
      double precision in_put(*), firfil(*), output(*)
      double precision alpha, beta
      integer incinp, i0_inp, n_inp
      integer incfir, i0_fir, n_fir
      integer incout, i0_out, n_out
c-----------------------------------------------------------------------------
c
c   dconv1d  performs a 1D convolution in the time domain :
c	O(j) = Beta*O(j) + Alpha * Sum[ I(i) * F(j-i) ]
c
c-----------------------------------------------------------------------------
c
c   PARAMETERS:
c
c	in_put:	Pointer to FIRST sample of sequence "in_put"
c	incinp:	Increment between two successive values of "in_put"
c	i0_inp: Index of the first element of "in_put" 
c	n_inp:	Number of samples of "in_put"
c
c	firfil:	Pointer to FIRST sample of sequence "firfil"
c	incfir:	Increment between two successive values of "firfil"
c	i0_fir: Index of the first element of "firfil" 
c	i0_fir:	Number of samples of "firfil"
c
c	output:	Pointer to FIRST sample of sequence "output"
c	incout:	Increment between two successive values of "output"
c	i0_out: Index of the first element of "output" 
c	n_out:	Number of samples of "output"
c
c	beta:	Scaling factor for the Output Sequence on Entry 
c	alpha:	Scaling factor for the convolution
c-----------------------------------------------------------------------------
c
c   PLEASE NOTE:
c	Please note that the array pointers must all point to the first 
c	element of the array "i0_inp", "i0_fir" and "i0_out".
c	 If "in_put" for example is defined as
c	    dimension in_put(-25:45)
c	Then "dconv1d" must be called with the following parameters
c	    call dconv1d( in_put(-25),1,-25,45, ... )
c
c-----------------------------------------------------------------------------
c
c   USAGE:
c	This module compute the result of the convolution in the "output" range
c	padding with zeroes when needed.
c
c	With some precautions it is possible that 
c		the Output sequence   OVERWRITE   the Input one.
c	Then (e.g. i0_out == i0_inp), it is necessary that   i0_fil >= 0
c	Said in Digital Signal Processing (DSP) jargon: "firfil" must be CAUSAL
c
c	In theory, an input sequence of "n_inp" samples starting at time "i0_inp", 
c	filtered by a sequence of "n_fir" samples starting at time "i0_fir",
c	will result in a new signal of (n_inp + n_fir - 1) non zero samples starting at
c	time (i0_inp + i0_fir). We just compute here the values that fall in that range
c	and zero the rest.
c	This may be interesting, for example when filtering a sequence of N samples, with a
c	symmetric filter of 2m+1 samples. If one wants only to compute the central N resulting
c	samples, the following call can be used:
c	    call dfir1d( f, 0, 1, N,  g, -m, 1, 2*m+1,  h, 0, 1, N)
c
c-----------------------------------------------------------------------------
c
c   This Fortran Subroutine written by
c	Jean-Pierre Panziera
c	Silicon Graphics Inc
c	September 20, 1991
c
c-----------------------------------------------------------------------------

      double precision tmp, zero, one
      parameter (zero = 0.0, one = 1.0)
      integer i1_inp, i1_fir, i1_out, k0, k1
      integer i, j, kout, kinp, kfir, n0, n1, jmin
c############################################################################-
c
c
c   Compute the FIR convolution
c
c
c############################################################################-
	i1_inp = i0_inp + n_inp - 1
	i1_fir = i0_fir + n_fir - 1
	i1_out = i0_out + n_out - 1
	k0 = i0_inp + i0_fir
	k1 = i1_inp + i1_fir

	n0 = i0_out
	if( n0 .lt. k0)		n0 = k0
	n1 = i1_out
	if( n1 .gt. k1)		n1 = k1
c-----------------------------------------------------------------------------
c
c	Scale the output
c
c-----------------------------------------------------------------------------
	kout = 1 
	do j = i0_out, i1_out
	    output(kout) = beta * output(kout)
	    kout = kout + incout
	end do
c-----------------------------------------------------------------------------
c
c	Compute the convolution
c
c-----------------------------------------------------------------------------
	kout = 1 + (n1 - i0_out) * incout
	do j = n1, n0, -1
	    tmp = zero
	    k0 = max(j- i1_fir, i0_inp)
	    k1 = min(j- i0_fir, i1_inp)
	    kinp = 1+ (k0-i0_inp) * incinp
	    kfir = 1+ (j-k0-i0_fir) * incfir
	    do i = k0, k1
		tmp = tmp + in_put(kinp) * firfil(kfir)
		kinp = kinp + incinp
		kfir = kfir - incfir
	    end do
	    output(kout) = output(kout) + alpha * tmp
	    kout = kout - incout
	end do

      return
      end
